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1.  INTRODUCTION 


This  report  covers  a  number  of  topics  related  to  the  application  of  finite  element  methods 
to  ballistic  structural  problems.  All  of  these  topics  are  extensions  and  modifications  of  finite 
element  formulations  and  computer  codes  previously  developed  in  this  technical  area. 

The  first  topic  to  be  discussed  involves  an  extension  of  so-called  gap  element  formulation 
(Zak  1 984).  In  a  previous  study,  a  finite  element  model  was  developed  which  can  handle 
material  contact  problems  inside  of  a  structure.  This  occurs  when  the  unloaded  structure 
contains  small  gaps  which  can  close  when  the  load  is  applied.  The  contact  of  the  adjacent 
portions  of  the  structure  leads  to  a  nonlinear  response.  The  finite  element  model  previously 
developed  for  this  problem  makes  use  of  so-called  gap  elements.  In  this  approach,  the  gaps 
are  modelled  by  these  elements  which  simulate  the  nonlinear,  physical  phenomena  of  material 
contact.  This  finite  element  model  was  incorporated  into  a  three-dimensional,  axisymmetric 
program  call  SAGA.  In  this  application,  the  gaps  were  allowed  to  exist  either  in  the  I  or  J 
directions  in  the  l-J  plane  of  the  finite  element  grid.  The  restriction  of  the  gap  geometry  to  be 
in  either  one  of  these  directions  placed  a  constraint  on  the  type  of  geometries  which  could  be 
handled.  For  example,  if  the  gap  made  a  number  of  sharp  changes  in  direction,  the  original 
formulation  woula  not  be  applicable.  Consequently,  it  was  decided  to  develop  a  new 
computer  code  for  the  contact  problem  in  which  this  restriction  would  be  removed.  The 
development  of  this  code  will  be  described  in  this  report. 

The  second  topic  of  this  report  will  deal  with  another  modification  to  the  original  gap 
element  formulation  (Zak  1984).  In  the  original  finite  element  formulation,  a  rather  simplified 
model  was  employed  to  represent  material  slip  after  gap  closure.  If  the  closure  of  the  gap 
occurs  during  the  loading  cycle,  the  gap  surfaces  will  not  move  relative  to  each  other  in  the 
direction  normal  to  the  gap.  However,  in  the  tangential  direction,  relative  motion  is  possible. 

In  the  original  gap  element  model,  the  tangential  forces  were  kept  constant.  The  objective  of 
the  modified  model  is  to  make  these  tangential  forces  proportional  to  the  normal  compressive 
force  across  the  gap.  This  feature  was  added  to  the  new  gap  element  code  which  allows  for 
combined  l-J  directions  of  the  gap. 
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The  third  topic  to  be  dealt  with  in  this  report  consists  of  adding  rate-dependent  plasticity  to 
the  gap  computer  code  (Drysdale  and  Zak  1983).  In  a  previous  study,  a  material  model  was 
developed  in  which  the  yield  condition  is  a  function  of  plastic  strain  rate.  This  model  was 
applied  to  finite  element  theory  and  incorporated  into  a  computer  program  called  SANX 
(Zak,  Craddock,  and  Drysdale  1979),  The  present  objective  was  io  use  this  rviaterial  model  m 
the  gap  element  computer  code.  In  addition,  some  further  improvements  to  this  model  were 
incorporated  so  as  to  eliminate  numerical  oscillations  which  occurred  in  SANX  when  abrupt 
changes  were  made  in  the  loading  conditions.  For  example,  one  interesting  situation 
consisted  of  a  slender  structural  member  subject  to  torsional  load  which  was  held  constant  at 
yield  and  then  followed  by  constant  normal  strain  rate.  This  loading  produces  large 
discontinuities  and,  consequently,  numerical  oscillations.  Reformulation  of  the  computational 
procedure  has  eliminated  numerical  oscillations. 

The  fourth  and  the  last  topic  to  be  discussed  involves  some  preliminary'  investigation  of 
code  combinations.  During  a  number  of  past  investigations,  many  different  but  related  finite 
element  codes  have  been  developed  for  ballistic  applications.  The  purpose  of  the  codes  is  to 
handle  a  number  of  different  configurations,  material  properties,  and  loading  conditions. 

During  the  present  investigation,  some  preliminary  work  was  done  on  the  feasibility  of 
combining  these  codes  so  that  advantage  can  be  taken  of  common  features  such  as  the  finite 
element  g'^'d  generation  and  other  input  preparation. 

2.  EXTENDED  GAP  GEOMETRY 

2.1  Introduction  to  SAGA  Program.  The  original  gap  element  computer  program  is  called 
SAAC  (Zak  1984).  The  starting  point  for  the  development  of  this  program  was  the  SAGA 
program.  Before  describing  the  SAAC  program  and  its  extension  in  this  investigation,  it  is 
useful  to  briefly  review  the  SAGA  code  and  its  application. 

Originally,  the  SAGA  program  was  developed  for  elastic  materials.  The  elastic  SAGA  itself 
is  an  extension  of  an  axisymmetric  code  called  SAAS.  The  SAAS  is  a  basic  finite  element 
program  with  two  degrees  of  freedom  at  each  node.  Some  time  ago,  it  was  useful  to  extend 
this  code  to  an  axisymmetric  code  with  three  degrees  of  freedom  at  each  node.  This  new 
code  was  denoted  as  SAGA.  The  motivation  for  this  new  code  was  to  permit  solution  of 
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problems  containing  torsional  loads  and  problems  involving  orthotropic  materials  where  the 
a-.ts  of  orthotropy  is  out  of  the  meridian  plane. 

Before  modifying  the  SAGA  program  to  incorporate  the  gap  element  model,  it  was 
necessary  to  develop  an  elastic-plastic  version  called  SAGAEP.  This  program  handles  the 
problem  by  incremental  solution.  Both  the  loads  and  the  materia!  properties  are  varied 
incrementally.  The  incremental  nature  of  the  SAGAEP  program  was  necessary  since  the  gap 
program  SAAC  is  nonlinear  even  if  the  material  properties  are  elastic.  The  nonlinear  nature 
arises  from  the  gap  contact  phenomena.  However,  one  of  the  original  requirements  in  the 
SAAC  development  was  to  incorporate  plasticihy. 

2.2  introduction  to  SAAC  Program.  Since  the  original  SAAC  program  forms  the  basis  for 
three  aspects  of  this  investigation,  it  is  necessary  to  describe  briefly  the  main  features  of  the 
gap  element  model  and  its  incorporation  into  the  computer  code. 

Tne  basic  idea  of  gap  elements  is  that  they  simulate  the  gap  before  closure  and  after 
Closure.  The  gap  element  material  properties  are  altered  to  simulate  material  contact.  This  is 
done  by  originally  assigning  very  low  stiffness  properties  to  the  element.  If  contact  occurs, 
certain  of  the  stiffness  properties  are  changed  to  make  the  element  very  stiff  in  the  direction 
perpendicular  to  the  gap  direction. 

The  gap  finite  element  model  can  be  described  by  referring  to  Figure  1  through  8. 

Figure  9  shows  the  flow  chart  for  the  program  SAAC.  The  basic  ideas  of  the  gap  element 
model  are  first  developed  with  reference  to  the  l-J  grid.  The  l-J  coordinate  system  will  locate 
elements  and  nodal  points.  To  illustrate  how  nodal  points  are  located,  take  a  typical  element 
in  an  l-J  grid  (Figure  1).  The  element  has  nodal  point  IX(N,1)  located  at  the  corner  that  lies 
closest  to  the  origin.  The  remaining  three  nodal  points  IX(N,2),  iX(N,3),  and  IX(N,4)  are 
obtained  by  moving  in  the  counterclockwise  direction  around  the  element.  The  direction  of  the 
gap  depends  on  an  input  parameter  called  the  principal  gap  direction.  As  an  example,  if  the 
principal  gap  direction  parameter  is  chosen  to  be  equal  to  1;  then,  the  gap  direction  will  be 
parallel  to  the  J-axis  while  the  direction  of  closure  will  be  parallel  to  the  l-axis  (Figure  2).  If 
the  principal  gap  direction  is  set  to  2,  the  direction  of  the  gap  will  be  parallel  to  the  l-axis  while 
the  gap  closure  will  be  in  the  J-axis  direction  (Figure  3).  In  the  present  model  and  the 
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program  SAAC,  there  is  only  one  choice  available  for  the  gap  direction  parameter.  For 
example,  if  the  parameter  is  chosen  to  be  1,  all  gap  elements  will  have  that  direction. 

The  gap  direction  will  determine  which  node  points  will  be  opposite  one  another  along  the 
gap.  This  information  is  important  for  reasons  that  will  be  best  understood  from  the  following 
discussion.  First,  the  program  must  know  which  node  points  to  pair  and  then  test  for  contact; 
and  secondly,  when  contact  is  established,  the  nodal  point  pairs  will  be  used  to  properly 
modify  gap  element  material  properties.  Referring  to  Figure  2,  if  the  principal  gap  direction  is 
set  to  1,  then  the  two  nodal  pairs  for  the  element  are  (IX[N,1],  IX[N,2])  and  (1X[N,3],  IX[N,4]). 
Similarly,  for  a  principal  gap  direction  of  2.  the  two  pairs  are  (IX[N,1],  IX[N,4])  and  (1X[N,2]. 
IX[N,3]), 


The  third  geometric  parameter  has  to  do  with  the  orientation  of  the  gap  element  relative  to 
the  R-2  axes.  This  orientation  will  be  used  to  transform  nodal  positions  and  nodal 
displacements  into  a  coordinate  system  where  one  axis  coincides  with  the  direction  of  the 
gap,  while  the  other  axis  coincides  with  the  direction  of  the  closure.  The  second  use  of  the 
gap  orientation  will  insure  that  gap  material  properties  will  only  stiffen  in  the  direction  of 
closure  once  con‘act  is  established. 

Once  the  principal  gap  direction  is  set,  the  angle  of  the  gap  can  be  calculated.  This  angle 
IS  defined  in  a  manner  indicated  by  Figure  4.  As  an  example,  take  a  typical  element  with  a 
principal  gap  direction  of  1  and  nodal  points  11,  12,  13,  and  14  [equal  to  iX(N,1),  IX(N,2), 
IX(N,3),  and  IX(N,4)  respectively].  The  angle  p  will  be  calculated  by  first  producing  a  "mean” 
line  through  the  element  (Figure  5).  The  coordinates  of  points  A  and  B  are,  respectively,  as 
follows: 


„  7  _  R{I1)  .  R{I4)  Z{I1)  .  Z{/4) 

A'  A  2  *  2 

R  ,Z  ^  ^(^2)  .  R(I3)  _  Z(/2)  ^  Z(/3) 
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Using  the  expressions  (Equations  1  and  2),  the  angle  p  can  be  calculated  by  using  the 
following; 


P  =  ARCtan 


(«B- 

(^B  -  Z,)  ■ 


(3) 


This  definition  of  p  will  have  coordinate  axis  N  coincident  with  the  direction  of  closure  and 
coordinate  axis  T  coincident  with  the  gap  direction. 

If  the  principal  gap  direction  is  equal  to  2,  the  angles  (Equations  1  and  2)  are  changed  by 
pairing  off  points  11  with  12  and  12  with  14.  Otherwise,  the  calculation  for  p  is  the  same  as  in 
Equatic.i  4.  The  establishment  of  the  nodal  pairs  also  permits  the  calculation  of  the  gap  size 
at  each  end  of  the  element.  The  significance  of  the  gap  size  is  that  it  is  later  used  to  check  if 
gap  closure  occurs.  In  the  SAAC  program,  the  geometrical  properties  of  gap,  including  the 
orientation  and  the  gap  size,  are  calculated  in  the  Subroutine  GAPANG  shown  in  Figure  9. 

The  solution  of  the  structural  problem  is  performed  incrementally.  The  external  load  is 
divided  into  load  steps.  At  each  load  step,  the  solution  checks  for  closure  at  each  set  of  nodal 
pairs  for  each  gap  element.  The  closure  is  determined  by  comparing  the  normal  displacement 
components  with  the  original  gap  size.  Determining  whether  or  not  opposing  node  points 
have  contacted  is  done  by  first  transforming  nodal  displacements  and  nodal  positions  into  the 
N,  T,  S  coordinate  system.  The  coordinates  N  and  T  are  given  in  Figure  4,  and  S  is  defined 
to  be  orthogonal  to  N  and  T.  Transformation  from  the  R,  Z,  0  coordinate  system  to  the  N,  T, 

S  coordinate  system  is  of  the  following  form; 


Un 

Ur 

•  •  1  T"!  ■ 

Uz 

.^s  . 

(4) 


where  C/^,  U^,  and  Uq  are  the  nodal  displacements  in  the  R,  Z,  and  0  directions,  respectively. 
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The  transformation  matrix  [  7]  is  as  follows; 


[  T\  = 


S/A/p  COSP  0 
-COSp  S/A/p  0 
0  0  1.0 


(5) 


where  the  angle  p  is  defined  by  Equation  3.  As  can  be  seen  in  Equation  5,  the  coordinate 
axes  0  and  S  will  always  coincide  in  this  transformation. 


The  transformation  of  nodal  positions  is  similar  to  the  transformation  of  nodal 
displacements.  This  transformation  is  of  the  following  form: 


(6) 


where  and  are  the  nodal  positions  in  the  R  and  Z  directions,  respectively.  P^  and  Pj 
are  the  nodal  positions  in  the  N  and  T  directions,  respectively.  The  transformation  matrix 

[  7'  ]  is  as  follows: 


[  r] 


S/A/p  COSp 
-COSp  S/A/p 


(7) 


Since  we  are  dealing  with  axisymmetric  elements,  the  Ug  does  not  affect  gap  closure.  The 
second  step  in  the  test  for  contact  involves  updating  nodal  positions  in  the  N  direction.  The 
updating  occurs  at  the  end  of  every  load  step.  This  procedure  is  applied  to  each  nodal  pair 
for  each  gap  element.  The  new  relative  position  of  the  nodal  pairs  determines  the  gap 
closure.  When  the  gap  is  closed,  no  further  relative  motion  in  the  normal  gap  direction  will  be 
allowed.  When  real  gaps  close,  the  two  opposing  material  faces  will  never  overlap  one 
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another.  For  an  effective  simulation  of  real  gaps,  this  behavior  must  be  duplicated  in  the  finite 
element  model.  It  is  conceivable  that  a  nodal  pair  that  is  not  contacting  at  the  end  of  a  given 
load  step  can  be  overlapping  at  the  end  of  the  next  load  step.  This  can  occur  because  the 
program  only  checks  for  contact  at  discrete  load  intervals  and  not  continuously.  The  analysis 
handles  the  problem  of  overlapping  by  splitting  load  steps.  The  first  step  in  this  process  takes 
all  the  nodal  pairs  that  are  overlapping  at  the  end  of  a  given  load  step  and  calculates  a  factor 
for  each  one  that  depends  on  the  degree  of  overlap.  This  procedure  works  something  like 
this:  When  a  particular  load  step  causes  a  nodal  pair  closure  and  an  overlap,  the  load  step  is 
modified.  For  example,  suppose  a  load  step  causes  nodal  pair  relative  displacement  twice  as 
large  as  necessary  to  just  close  the  gap.  Then  this  load  step  is  modified  by  a  factor  of 
one-half,  and  the  analysis  is  repeated  for  the  new,  smaller  load.  Testing  for  contact  and  the 
calculation  of  load  factors  occur  in  Subroutine  TEST.  The  Subroutine  TEST  is  called  from 
Subroutine  SOLV  after  nodal  displacements  have  been  calculated.  An  assumption  is  made 
to  the  load  factor  analysis  to  reduce  execution  time.  The  assumption  is  that  the  contact 
condition  of  nodal  pairs  is  determined  by  the  original  load  step  and  not  by  the  reduced  load 
step.  As  an  example  of  this,  if  a  given  number  of  nodal  pairs  are  overlapping  at  the  end  of  a 
load  step,  the  program  will  then  assume  that  these  nodal  pairs  will  remain  in  contact  even 
though  the  load  step  will  be  reduced.  An  advantage  of  this  assumption  is  that  original  load 
steps  are  only  split  once.  A  disadvantage,  however,  is  that  the  program  might  predict  gap 
closure  at  a  different  load  as  compared  to  a  real  structure  with  the  same  configuration.  This 
problem  is  minimized  by  keeping  load  steps  small. 

After  gap  closure  has  been  established,  the  nodal  pairs  involved  will  not  move  relative  to 
each  other  in  the  N  direction.  Unloading,  of  course,  could  separate  the  gap,  and  program 
SAAC  can  handle  this  possibility.  The  gap  element  simulates  the  contact  situation  by 
changing  certain  stiffness  properties  of  the  element.  This  is  not  done  by  changing  the 
material  properties  but  rather  by  directly  changing  the  stiffness  matrix  of  the  gap  element. 

The  global  coordinates  which  the  SAAC  program  uses  are  the  cylindrical  coordinates  R,  Z, 
and  0.  Consequently,  the  stiffness  matrix  for  all  elements  is  assembled  in  these  coordinates. 
However,  in  order  to  simulate  the  gap  closure,  it  is  necessary  to  change  the  element  stiffness 
in  the  N,  T,  S  coordinate  system.  The  element  stiffness  matrices  in  the  SAGA  and  the  SAAC 


7 


programs  represent  12  degrees  of  freedom.  Suppose  the  stiffness  matrix  in  the  N,  T,  S 
system  is  represented  by  the  following: 

[S]  =  .  (8) 

/=  1.12 

Before  contact,  the  material  properties  of  the  element  are  chosen  to  be  very  soft  relative  to 
the  structural  material  surrounding  the  gap.  When  the  stiffness  matrix  (Equation  8)  is 
calculated  for  the  gap  element,  the  stiffness  coefficients  S,y  will  be  small  relative  to  the 
adjacent  structure.  For  numerical  reasons,  the  material  properties  cannot  be  identically  zero. 
When  contact  is  established,  the  stiffness  coefficients  in  the  N,  T,  S  system  are  modified  to 
simulate  the  fact  that  two  of  the  nodal  points  will  move  together.  This  is  achieved  by 
increasing  the  numerical  size  of  certain  of  the  stiffness  coefficients.  Theoretically,  the 
increasing  value  of  these  coefficients  could  be  infinite;  however,  it  has  been  found  that  this 
leads  to  numerical  difficulties.  Consequently,  the  program  SAAC  increases  these  coefficients 
so  that  they  are  approximately  two  magnitudes  larger  than  similar  coefficients  for  adjacent 
structural  materials. 

The  procedure  of  changing  the  stiffness  and  simulating  contact  across  the  gap  can  be 
illustrated  as  follows. 

A  gap  element  has  four  different  contact  states  for  each  value  of  the  principal  gap 
direction.  Consider  an  element  with  a  principal  gap  direction  of  1  and  opposing  nodal  pairs 
(1,2)  and  (3,4)  (Figure  7).  The  following  four  contact  states  are  possible: 

(1)  State  No.  1  -  both  node  pairs  not  contacting 

(2)  State  No.  2  -  (1,2)  contacting;  (3,4)  not  contacting 

(3)  State  No.  3  -  (1,2)  not  contacting;  (3,4)  contacting 

(4)  State  No.  4  -  both  nodal  pairs  contacting. 

The  gap  element  stiffness  matrix  is  modified  differently  for  each  of  the  four  states.  The 
proper  modification  of  the  element  stiffness  matrix  is  as  follows: 
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(1)  State  No.  1  -  all  matrix  elements  remain  unchanged  and  small 

(2)  State  No.  2  -  set  S, ,  =  =  -S,  ^  =  -S^  ,  =  large  number 

(3)  State  No.  3  -  set  S77  =  =  -S,o,7=  large  number 

(4)  State  No.  4  -  set  Sj  ^  ^  y  =  Siq  ^q  =  -S«  =  ,  =  'Sy  =  -S  ^qj 

=  large  nurhber. 

If  the  principal  gap  direction  of  the  element  is  equal  to  2,  the  gap  element  stiffness  matrix 
modification  will  be  similar  to  that  shown  previously.  However,  different  stiffness  coefficients 
will  be  changed  as  indicated  by  appropriate  degrees  of  freedom  for  the  nodal  pairs. 

After  modification  of  the  stiffness  coefficients  in  the  N,  T,  S  coordinate  system,  the 
stiffness  matrix  is  transformed  into  the  cylindrical  coordinates.  The  matrix  transformation 
procedure  is  given  by  the  following: 

[K]  =  [TV  [S]  [T]  ,  (9) 

where  [S]  and  [K]  are  the  element  stiffness  matrices  in  the  N,  T,  S  and  R,  Z,  0  coordinate 
systems,  respectively.  The  transformation  matrix  [T]  is  as  follows: 

[t] 

[t] 

[t] 

[t] 


where  [f]  is  given  by  the  following: 

s//vp  -COSp  0 

[f]  =  COSp  S/A/p  0  (11) 

0  0  1.0 


9 


In  the  computer  program  SAAC,  the  modification  of  the  gap  element  stiffness  matrix  and 
the  coordinate  transformation  are  performed  in  Subroutine  RETRANS  as  shown  in  Figure  9. 

2.3  Extended  General  Gap  Model.  The  program  SAAC,  with  the  gap  element  model 
described  previously,  forms  the  starting  point  for  a  modified  program  called  SAAC2.  This 
program  contains  the  same  basic  ideas  as  SAAC  but  allows  for  a  general  gap  configuration. 
This  configuration  permits  the  gap  orientations  to  exist  in  I  and  J  directions  simultaneously. 

To  illustrate  this,  consider  a  corner  gap  in  the  1  and  J  coordinates  as  shown  in  Figure  10. 

This  figure  shows  three  gap  elements.  Element  No.  1  is  in  the  I  direction;  element  No.  2.  is  in 
the  J  direction  and  the  gap  element.  No.  3,  is  a  corner  gap  element.  For  these  gaps  for  which 
this  model  is  designed,  the  corner  element  represents  a  very  small  volume  of  the  structure 
and,  therefore,  is  of  negligible  contribution.  Consequently,  in  the  analysis  of  SAAC2  corner 
gap  elements  such  as  No.  3  in  Figure  10,  always  retain  constant,  low  stiffness. 

The  program  SAAC2  retains  the  same  arrangement  as  SAAC  and,  therefore,  has  the 
same  flow  chart  as  shown  in  Figure  9.  The  basic  change  in  the  new  model  is  to  allow  each 
element  to  possess  its  own  direction  parameter.  This  is  achieved  by  adding  this  parameter  to 
material  block  input  cards.  For  normal  elements,  not  gap,  this  parameter  is  left  blank  on  the 
input  and  is  ignored. 

In  order  to  create  program  SAAC2  from  SAAC,  changes  have  been  made  in  the  following 
parts  of  the  program:  MAIN,  GAPANG,  POINTS,  SOLVE,  and  STIFF. 

In  Subroutine  POINTS,  one  material  block  parameter  IJDIR(M)  is  read  in  for  each  block. 
Using  this  block  parameter,  an  element  parameter  IG(N)  is  set  up  for  each  element  in  the 
structure.  As  mentioned  before,  this  parameter  is  automatically  taken  equal  to  zero  if  no  input 
in  the  material  block  is  entered.  This  applies  to  nongap  elements.  In  the  original  SAAC 
program,  there  was  only  one  gap  direction  parameter— LGAP.  This  parameter  has  been 
retained  in  SAAC2  problem;  however,  it  is  continuously  reset  for  each  individual  element. 

The  input  instructions  for  SAAC2  are  described  in  Appendix  A.  These  instructions  can 
also  be  used  for  SAAC  since  the  differences  are  small;  and  when  these  differences  occur, 
they  are  indicated. 
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The  program  SAAC2  was  checked  out  on  a  number  of  test  examples.  One  of  the  initial 
checks  performed  was  to  compare  numerical  answers  with  those  from  SAAC  for  one- 
directional  gaps.  A  number  of  examples  for  SAAC  were  described  in  the  previous  report  (Zak 
1984). 

3.  GAP  MODEL  INCLUDING  SLIP 


3.1  Background  and  Objectives.  In  the  original  study  which  resulted  in  the  development 
of  the  SAAC  program  (Zak  1984),  a  preliminary  attempt  was  made  to  account  for  material  slip 
forces  after  the  gap  was  closed.  As  a  result  of  this  effort,  a  separate  version  of  the  program, 
called  SAACFR,  was  developed  which  attempted  to  include  slip  forces  after  material  contact 
was  achieved.  This  model  was  rather  simple  and  was  designed  only  to  establish  the 
feasibility  of  including  such  effects  in  the  gap  element  model  One  of  the  objectives  of  the 
present  study  was  to  develop  a  more  comprehensive  slip  model  which  would  be  more 
realistic.  Before  describing  the  new  slip  model,  it  is  useful  to  review  briefly  the  slip  model  in 
the  original  SAACFR  program. 


In  the  original  SAAC  program,  no  allowance  was  made  for  the  tangential  forces  due  to 
material  slip  after  gap  closure  was  achieved.  In  order  to  investigate  ♦he  relative  effect  of  this 
phenomena,  a  simple  slip  model  was  introduced.  The  main  features  of  this  model  are  as 
follows.  The  model  introduces  frictional  forces  which  modify  the  gap  stiffness  model.  The 
effects  of  friction  are  related  to  the  stiffness  in  the  direction  of  closure.  This  relationship  is 
made  through  the  use  of  a  frictional  coefficient.  As  an  example,  if  nodal  pair  (1,2)  in  an 
element  is  contacting,  then  the  following  relationships  will  model  friction: 


^2.2  -  ^5,5  -  '^5.2  -  '^2.5  “  V 
^3.3  ~  ^6.6  ~  '^3,6  ~  '^6,3  ~ 


(12) 

(13) 


where  ^  is  a  frictional  coefficient  or  simply  a  proportionality  constant  between  the  stiffness 
coefficients  perpendicular  to  the  gap  and  in  the  tangential  directions.  The  two  normal  degrees 
of  freedom  are  1  and  4  and,  therefore,  the  four  tangential  degrees  of  freedom  are  2,  3,  and 
5,  6.  By  examining  Equations  12  and  13,  it  can  be  seen  that  they  relate  the  tangential 
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stiffness  to  the  normal  stiffness  after  gap  closure.  In  that  sense,  this  modification  does  restrict 
the  tangential  motion;  however,  it  does  not  fully  represent  what  is  normally  considered  as  a 
frictional  behavior.  In  spite  of  the  limitations  of  this  model,  it  was  useful  to  examine  the 
relative  effect  of  the  slip  on  the  response  of  a  number  of  test  examples.  Equations  12  and  13 
can  be  extended  easily  to  nodal  pair  (3,4).  If  the  nodal  pairs  are  (2,3)  and  (1,4),  similar 
modified  relations  can  be  obtained.  The  SAACFR  program,  since  it  is  an  extension  of  SAAC 
code,  automatically  determines  the  nodal  pairs  and  applies  the  tangential  correction  as 
illustrated  by  Equations  12  and  13.  In  the  test  examples  executed,  the  value  of  frictional 
coefficient  ^  had  to  be  estimated  based  on  some  reasonable,  preliminary  estimates. 

The  modification  of  stiffness  and  the  friction  analysis  are  both  done  in  Subroutine 
RETRANS.  This  change  of  gap  material  properties  takes  place  in  the  N,  T,  S  coordinate 
system  and  afterwards  is  transformed  back  in  the  R,  Z,  0  system. 

3.2  Extended  Slip  Model.  As  pointed  out,  the  original  slip  model  was  designed  to  be 
preliminary  and  did  not  possess  properties  which  normally  are  associated  with  frictional 
phenomenon  of  sliding  surfaces.  The  new  mode!  which  has  been  developed  attempts  to 
overcome  these  objectives. 

The  force  in  slip  phenomena  is  represented  by  Colomb  friction.  This  can  be  illustrated  by 
force-displacement  relation  as  shown  in  Figure  11.  Up  to  a  certain  point,  the  force  and 
displacement  are  linearly  related  and  after  certain  value  of  force  is  exceeded,  the  force  will 
remain  reasonably  constant.  It  can  be  seen  from  Figure  1 1  that  the  forces  involved  in  the  slip 
phenomena  have  approximately  the  same  behavior  as  an  elastic  perfectly  plastic  stress-strain 
relation.  In  an  elastic  perfectly  plastic  response,  the  stress  at  which  the  stress-strain  curve 
becomes  flat  is  known  as  the  yield  stress.  In  the  case  of  the  frictional  phenomena,  the  yield 
stress  would  correspond  to  the  force  at  which  the  force  displacement  becomes  flat.  However, 
unlike  in  tne  yield  phenomena  where  yield  stress  is  relatively  constant,  the  value  of  force 
where  the  curve  becomes  flat  (Figure  11)  depends  on  the  normal  forces.  Consequently,  the 
slip  phenomena  can  be  modeled  on  a  basis  of  an  analogy  with  elastic-plastic  response  where 
the  yield  stress  becomes  a  function  of  normal  loads  to  the  direction  of  the  gap.  This  approach 
has  been  taken  in  the  development  of  an  improved  slip  model. 
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The  modified  slip  model  involves  the  modification  of  certain  of  the  stiffness  coefficients  in 
the  gap  element.  Consider  now  the  modification  of  Equations  12  and  13  for  nodal  pair  (1,2). 

S2  2  =  S5  5  =  -i>52  =  '^2,5  ~ 

^3,3  ~  ^6.6  -  ~^3,6  ~  '^6,3  ~  ^  5) 

where  in  Equations  14  and  15,  S  represents  the  value  of  the  stiffness  coefficients  in  tangential 
directions  after  closure.  The  behavior  of  S  is  assumed  to  represent,  symbolically,  the  force- 
displacement  properties  as  expressed  by  Figure  11.  The  stiffness  coefficients  in  Equations  14 
and  1 5  represent  the  incremental  stiffness  relations.  The  S  parameter  is  related  to  the 
normal  forces  in  the  nodal  pairs  (1,2).  These  nodal  forces  have  the  degrees  of  freedom 
I  and  4,  respectively.  These  forces  are  equal  and  opposite.  The  actual  calculation  of  the 
parameter  S  and  its  relation  to  the  stresses  follows  the  following  procedure.  Consider  the 
nodal  pairs  (1 ,2).  Not  including  the  gap  element,  these  nodal  pairs  will  be  connected  to  four 
real  elements — two  elements  for  each  pairing.  The  special  case  is,  of  course,  when  nodal 
pairs  are  on  the  boundary,  in  which  case,  there  will  be  only  one  element  associated  with  each 
node.  The  modifications  expressed  by  Equations  14  and  15  take  affect  after  contact  between 
nodal  pairs  has  been  established.  During  the  first  load  steps  after  the  contact  has  been 
established,  the  parameter  S  is  chosen  to  be  large  and  of  the  same  magnitude  as  the 
modifications  to  stiffness  coefficients  corresponding  to  the  normal  directions.  This  procedure 
established  that  originally,  after  contact,  no  slip  occurs,  in  the  subsequent  load  steps,  the 
fallowing  check  for  slip  is  followed.  Using  the  elements  in  question,  the  computer  program 
evaluates  both  normal  stresses  and  the  shear  stresses  in  the  N,  T,  S  coordinate  system. 

These  stresses  are  then  averaged  over  the  four  or  two  elements.  Using  the  two  averaged 
shear  stresses,  a  resultant  shear  stress  is  calculated.  The  relation  between  the  average 
normal  stress  and  the  average  resultant  shear  stress  is  now  used  for  checking  and 
establishing  slip  conditions.  The  ratio  between  the  average  resultant  shear  and  the  normal 
stress  is  calculated.  If  this  ratio  is  less  than  input  value  of  friction  coefficient,  no  slip  is 
assumed  to  occur,  and  the  parameter  S  is  retained  at  a  large  value.  If  the  stress  ratio 
exceeds  the  prescribed  friction  coefficient,  the  slip  will  occur,  and  this  will  be  achieved  by 
reducing  the  stiffness  parameter  S  to  a  low  value.  The  low  level  of  S  is  used  to  simulate  the 
flat  response  in  Figure  11.  In  the  calculation,  it  was  found  that  the  low  value  of  S  could  be 
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chosen  to  be  of  order  of  magnitude  corresponding  to  the  original  gap  element  properties 
before  contact.  Similar  procedure  applies  for  nodal  pairs  (3,4).  For  the  second  orientation  of 
the  gap  in  the  l-J  coordinate  system,  the  nodal  pairs  involved  are  (1,4)  and  (2,3).  Similar 
relations  are  automatically  calculated  for  this  combination. 

The  procedure  for  the  slip  model  allows  for  additional  features.  For  example,  it  is  possible 
to  make  the  friction  coefficient  to  be  a  function  of  the  stresses  in  the  N,  T,  S  coordinate 
system.  Also,  as  the  slip  progresses,  it  is  possible  to  modify  the  first  part  of  the  curve  of 
Figure  11  to  be  stress  or  deformation  related.  The  deformation,  in  this  case,  would  be 
specified  by  amount  of  slip  that  has  already  occurred. 

4.  RATE  DEPENDENT  PLASTICITY 

4.1  Background  and  Objectives.  One  of  the  additional  objectives  of  the  present 
investigation  was  to  introduce  the  theory  of  rate-dependent  plasticity  into  the  SAAC  program 
(Drysdale  and  Zak  1983).  Since  the  present  investigation  has  developed  a  more  general  gap 
element  program,  SAAC2,  the  rate-dependent  model  was  introduced  into  the  SAAC2  code. 

The  rate-dependent  model  was  first  incorporated  into  the  finite  eletnent  formulation  in 
connection  with  the  SANX  computer  code  (Zak,  Craddock,  and  Drysdale  1979).  The  SANK 
code  represents  a  family  of  related  codes  which  has  been  developed  for  the  purpose  of 
approximate,  three-dimensional  finite  element  analysis  of  certain  configurations.  The 
configurations  for  which  SANX  programs  were  designed  involve  situations  which  have  an 
axisymmetric  appearance  but  have  nonaxisymmetric  features  which  render  a  purely 
axisymmetric  solution  inapplicable.  A  series  of  SANX  programs  was  developed  which 
included  elastic,  elastic-plastic,  elastic-viscoplastic,  and  rate-dependent  plastic  material 
versions. 

Before  describing  the  application  of  the  rate-dependent  plasticity  to  the  SAAC  program,  it 
is  useful  to  review  this  material  model  and  its  application  to  incremental  analysis.  In  its  simple 
form,  the  rate-dependent  plastic  model  can  be  described  by  the  following  yield  function: 

F(o,,e^,e'’)  =  f(S,  -  a,)  -  /<(e'>)  =  0  ,  (16) 
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where 


f(Sij  -  a,^  =  (Sjj  -  aij)  (S,j  -  a^j) 

Sjj  =  deviatoric  stress 

K=iA[a^{enf 

Oy  =  rate-dependent  yield  stress 

Ep  =  y%e/ye^  ,  which  is  defined  as  the  effective  plastic  strain  rate 
=  plastic  strain 

da,y  =  cdejj,  where  c  is  the  constant  kinematic  strain  hardening  parameter. 


The  empirical  relation  used  to  define  the  rate-dependent  yield  stress  is  as  follows: 


" 

f  ^ 

” 

1  +  b/n 

1 

_ 

_ 

(17) 


where 

Oq  =  static  yield  stress 

Eg  =  transition  strain  rate 

b  =  strain  rate  hardening  parameter. 


Using  Equations  16  and  17,  it  is  possible  to  develop  a  set  of  elastic-plastic,  incremental 
stress-strain  relations.  These  incremental  stress-strain  relations  are  then  used  in  the  finite 
element  model.  The  details  of  obtaining  the  incremental  stress-strain  relations  can  be  found 
elsewhere  (Drysdale  and  Zak  1 983)  and  will  not  be  repeated  here.  However,  it  is  useful  to 
summarize  the  steps  which  lead  from  Equations  16  and  17  to  the  incremental  relations. 
These  main  steps  are  as  follows: 


Step  1  During  the  yield  process,  the  condition  is  imposed  that  the  incremental 

changes  satisfy  the  requirement  that  the  yield  function  F  has  zero  change. 
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Step  2  Incremental  changes  in  plastic  strain  are  given  by  the  associated  plastic  flow 
rule. 

Step  3  Steps  1  and  2  above  yield  a  differential  equation  for  the  effective  plastic  strain. 

Step  4  For  each  time  increment,  a  solution  of  the  differential  equation  in  Step  3  is 
performed  to  obtain  the  effective  plastic  strain. 

Step  5  Step  4  permits  the  calculation  of  the  change  of  the  effective  plastic  strain  and 
the  change  in  the  actual  plastic  strain  during  each  interval. 

Step  6  Using  the  results  of  Step  5,  it  is  then  possible  to  write  the  incremental  stress- 
strain  relations. 

The  incremental  stress-strain  relations  are  obtained  by  applying  Steps  1  to  6.  The  general 
form  of  these  relations  is  as  follows: 


>  e-^-)  .  (18) 

def 

where  the  parameters  D.  A,  and  X  are  defined  during  the  development  of 

Equation  18.  It  may  be  noted  that  the  right-hand  side  of  Equation  18  forms  the  basis  for 
converting  the  SAAC2  program  to  include  the  rate-dependent  yield  stress  into  the  gap  element 
program. 

4.2  Modification  of  the  SAAC  Program.  The  rate-dependent  plastic  model  was 
incorporated  into  the  SAAC2  code.  The  new  version  of  the  code  is  called  SAAC3.  The 
changes  to  the  SAAC2  program  involved  the  following  parts  of  the  program:  MAIN,  ELPLSS, 
STIFF.  STRESS,  and  TRISTF. 

To  create  SAAC3  from  SAAC2  code,  it  is  first  necessary  to  incorporate  rate  calculations. 
The  rate  calculations  involve  the  evaluation  of  stress  rates  for  each  element  as  well  as  the 
rate  of  change  of  the  effective  plastic  strain  rate.  These  rates  are  calculated  in  SAAC3  by 
evaluating  the  changes  of  these  quantities  at  each  time  or  load  incremental  and  then  updating 
the  total  rates.  The  stress  rates  are  needed  in  the  calculation  of  the  parameter  A  in 
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Equation  1 8.  The  rate  of  the  effective  plastic  strain  rate  affects  the  derivative  of  K  relative  to 
The  bulk  of  the  calculations  involved  in  the  rate  calculations  are  included  in  the 
Subroutine  STRESS. 

The  second  major  change  involves  calculation,  at  each  interval,  of  extra  body  forces. 
From  Equation  18,  it  can  be  seen  that  the  second  term  on  the  right-hand  side  does  not 
contain  either  total  incremental  stress  or  strain.  Consequently,  when  Equation  18  is  used  in 
the  finite  element  model,  the  effect  of  this  term  is  to  produce  terms  which  act  as  body  forces. 
The  contribution  of  these  effective  body  forces  is  added  to  the  incremental  force  matrices. 
This  step  is  performed  in  the  Subroutine  TRISTF. 

The  actual  evaluation  of  the  second  term  on  the  right  side  of  Equation  18  involves  a 
number  of  different  calculations.  These  calculations  depend  on  the  stress  and  strain  rates, 
and  they  are  performed  in  the  Subroutine  ELPLSS.  In  order  to  create  SAAC3  from  the 
SAAC2  code,  this  subroutine  had  to  be  completely  rewritten.  The  basis  for  the  new 
Subroutine  ELPLSS  was  the  SANX  program  with  the  rate-dependent  yield  feature.  Starting 
with  SANX,  the  Subroutine  ELPLSS  was  modified  and  a  number  of  changes  made  in  the 
numerical  calculations.  These  modifications  corrected  some  numerical  stability  problems 
originally  observed  in  the  SANX  program.  The  new  version  of  the  Subroutine  ELPLSS  is 
given  in  Appendix  B. 

The  numerical  instability  originally  observed  relative  to  SANX  program  involved  special 
loading  situations.  The  particular  load  situation  to  which  the  SANX  code  was  applied,  as  a 
test  example,  involved  two  different  load  cycles  with  step  inputs.  First,  a  uniaxial  structural 
member  was  subject  to  a  torsional  load  until  yield  occurred.  After  achieving  yield,  the 
structure  was  subjected  to  a  constant  normal  strain  rate  load.  This  loading  subjects  the 
structure  to  a  step  load.  The  original  SANX  formulation  led  to  numerical  oscillations  as  the 
load  was  changed  from  torsion  to  uniaxial  load.  The  new  vefsion  of  the  Subroutine  EPLSS 
does  not  have  the  same  difficulty.  In  order  to  check  the  new  formulation  for  numerical 
instability,  a  special  version  of  SAAC3  was  prepared.  This  version  is  called  SAAC3X.  The 
reason  for  this  new  version  is  that  it  simulates  a  constant  strain  rate  loading.  Normally,  finite 
element  programs  allow  tor  either  force  input  or  displacement  boundary  conditions.  The 
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additional  version,  SAAC3X,  simulates  constant  strain  rate  through  displacement  boundary 
conditions. 

4.3  Input  Procedure.  The  various  gap  element  programs  developed  in  this  investigation 
have  similar  input  procedures.  The  description  of  the  input  data  for  all  these  programs  is 
similar  to  that  needed  for  the  original  SAAC  program.  Appendix  A  gives  the  input  procedure 
for  the  computer  programs  SAAC2  and  SAAC3.  It  may  be  noted  that  the  procedure  for 
SAAC2  and  SAAC3  is  similar  except  for  the  fact  that  SAAC3  has  a  new  input  card  in  the 
beginning.  This  card  is  labelled  as  Card  No.  0  in  Appendix  A.  There  is  no  difference  in  the 
input  procedures  for  SAAC3  and  SAAC3X.  However,  a  special  set  of  data  was  developed  for 
use  in  connection  with  the  special  loading  consisting  of  torsion  followed  by  normal  constant 
strain  rate.  This  data  was  used  to  check  numerical  stability  and  test  the  previous  results  from 
the  SANX  code.  The  same  results,  minus  the  numerical  instability,  were  achieved. 

5.  CODE  UNIFICATION 

5.1  Background.  The  last  topic  to  be  covered  in  this  report  involves  a  preliminary 
feasibility  study  involving  code  unification.  The  motivation  for  this  study  is  the  fact  that  a 
number  of  related  finite  element  codes  have  been  developed  with  similar  flow  charts  and 
common  parts.  These  programs  include  SAAS,  SAGA,  and  SANX  and  now  a  new  series  of 
SAAC  codes  described  in  this  report.  Since  these  codes  are  similar,  it  is  very  natural  to 
consider  if  they  can  be  combined  into  a  code  with  multiple  functions.  One  aspect  of  this 
study,  therefore,  was  involved  with  examining  this  question. 

5.2  Preliminary  Study.  A  quick  review  of  the  various  codes  makes  it  obvious  that  there 
are  some  easily  identified  common  aspects.  The  obvious  similarities  are  in  the  area  of  input 
procedure  and  in  the  fact  that  some  of  the  subroutines  are  basically  the  same  in  all  these 
codes.  The  main  similarity  lies  in  the  area  of  generating  the  finite  element  geometrical 
properties.  Consequently,  the  question  which  was  investigated  in  this  study  was  whether  the 
geometrical  portion  of  these  codes  could  be  separated  and  executed  as  a  separate  program. 
This  was  achieved  by  creating  a  truncated  version  of  a  finite  element  code  which  includes  the 
following  parts;  MAIN,  MESH,  POINTS,  MNIMX,  ANGLE,  and  CIRCLE.  This  new  program 
generates  the  finite  element  grid  properties.  Included  is  the  nodal  data,  element  information. 
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boundary  conditions,  and  surface  loads.  The  program  does  not  perform  any  solutions  since  it 
does  not  create  stiffness  properties,  matrix  assembly,  or  solution  steps.  In  that  sense,  the 
new  program  is  relatively  independent  of  any  of  the  original  codes.  One  application  of  this 
new  code  has  been  computer  plots  of  finite  element  geometries.  In  this  use,  the  code  has 
been  coupled  with  plotting  routines,  and  it  has  been  used  to  examine  finite  element  grids 
before  any  solution  is  attempted.  This  application  has  been  found  to  be  most  useful  in  grid 
generation. 

Future  plans  in  this  area  involve  generating  SAAS,  SAGA,  SANX,  and  SAAC  codes 
without  the  grid  generation  capability  and  driving  these  with  a  common  code  similar  to  the  one 
which  has  been  developed. 
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Gap  Direction 


Figure  5.  Gap  Element  Mean  Line  Used  to  Determine  Orientation. 
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Figure  6.  Degrees  of  Freedom  for  a  Typical  Element. 
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Figure  7.  Gap  in  the  J  Direction  Used  to  Illustrate  Various  Stiffness  Modifications. 


Figure  8.  Gap  in  the  I  Direction  Used  to  Illustrate  Various  Stiffness  Modifications. 
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Figure  9.  Flow  Chart  for  SAAC  Family  of  Codes. 
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Figure  11.  Force-Displacement  Relation  for  Material  Slip. 
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APPENDIX  A: 

INPUT  PROCEDURE  FOR  SAAC  PROGRAMS 
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DATA  NO.  0  -  RATE  PROPERTY  CARD 


(Needed  only  for  SAAC3  and  SAAC3X  Codes) 


Format  (5E12.6) 


Columns 


1-12  At  (Time  increment) 

13-24  b  (Rate  parameter  in  Equation  17) 

25-36  ^  (Rate  parameter  in  Equation  17) 

(Next  two  parameters  needed  only  for  SAAC3X) 
37-48  STRRATE  (Normal  strain  rate) 

49-60  STRSTEP  (Shear  stress  step  input) 


DATA  NO.  1  -  TITLE  CARD 
Format  (20A4) 

Columns  1-80  Title  (Title  for  particular  case) 
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DATA  NO.  2  -  CONTROL  CARD 
Format  (615.  F5.0.  615) 
Columns  1-5 


6-10 

11-15 

16-20 

21-25 

26-30 

31-35 

36-40 


41-45 

46-50 


51-55 


56-60 


NNLA  (Number  of  nonlinear 
approximations:  NNLA  =  1  for  this 
version  of  the  program.) 

NUMTC  (Number  of  temperature  cards; 
if  -2.  a  constant  temperature  is 
specified.) 

NUMMAT  (Number  of  different 
materials;  6  maximum) 

NUMPC  (Number  of  boundary  pressure 
cards;  200  maximum) 

NUMSC  (Number  of  boundary  shear 
cards;  200  maximum) 

NUMST  (Number  of  boundary  shear 
cards  in  tangential  direction; 

200  maximum) 

TREF  (Reference  temperature) 

INERT  (This  parameter  decides  if 
inertia  loads  will  be  present; 

INERT  =  0  means  zero  values  of 
axial  acceleration  and  angular 
acceleration  and  velocity  for 
each  load  increment.) 

NLINC  (Number  of  load  increments 
with  time;  NLINC  >  1) 

INCI  (If  INCI  =  0,  then  inertia 
loads  for  each  time  increment  will 
be  equal  to  the  inertia  load  for 
the  first  time  increment.) 

INCF  (If  INCF  =  0,  then  surface 
loads  for  each  time  increment  will 
be  the  same  as  for  the  first  time 
increment.) 

IPLOT  (Plot  parameter;  IPLOT  =  1 
if  plot  is  required.) 
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DATA  NO.  3  -  MESH  GENERATION  CONTROL  CARD 


Format  (515) 

Columns 

1-5 

MAXI  (Maximum  value  of  1  in  mesh; 

25  maximum) 

6-10 

MAXJ  (Maximum  value  of  J  in  mesh; 
100  maximum) 

11-15 

NSEG  (Number  of  line  segment  cards) 

16-20 

NBC  (Number  of  boundary  condition 
cards) 

21-25 

NMTL  (Number  of  material  block 
cards) 
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DATA  NO.  4  -  LINE  SEGMENT  CARDS 


The  order  of  line  segment  cards  is  immaterial,  except  when 
plots  are  requested;  in  this  case,  the  line  segment  cards 
must  define  the  perimeter  of  the  solid  continuously.  The 
order  of  line  segment  cards  defining  internal  straight  lines 
is  always  irrelevant. 


FORMAT  (3[213,  2F8.3],  15) 
Columns  1-3 

4-6 
7-14 
15-22 
23-25 
26-28 
29-36 
37-44 
45-47 
48-50 
51-58 
59-66 
67-71 


I  coordinate  of  the  1st  point 
J  coordinate  of  the  1st  point 
R  coordinate  of  the  1st  point 
Z  coordinate  of  the  1st  point 
I  coordinate  of  the  2nd  point 
J  coordinate  of  the  2nd  point 
R  coordinate  of  the  2nd  point 
Z  coordinate  of  the  2nd  point 
I  coordinate  of  the  3rd  point 
J  coordinate  of  the  3rd  point 
R  coordinate  of  the  3rd  point 
Z  coordinate  of  the  3rd  point 
Line  segment  type  parameter 


If  the  number  in  coiumn  71  is 

0  Point  (input  only  1st  point) 

1  Straight  line  (input  only  1st 
and  2nd  points) 

2  Straight  line  as  an  internal 
diagonal  (input  only  1st  and  2nd 
points) 
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Circular  arc  specified  by  1  st  and 
3rd  points  at  the  ends  of  the  arc 
and  2nd  point  at  the  midpoint  of 
the  arc 

4  Circular  arc  specified  by  1  st  and 
2nd  points  at  the  ends  of  the  arc 
with  the  coordinates  of  the  center 
of  the  arc  given  as  the  3rd  point 
(delete  I  and  J  for  3rd  point) 

5  Straight  line  as  boundary  diagonal 
for  which  I  of  1st  point  is  minimum 
for  its  row  and/or  I  of  2nd  point 

is  minimum  for  its  row  (input  only 
1  st  and  2nd  points) 

6  Straight  line  as  boundary  diagonal 
for  which  I  of  1st  point  and/or 
2nd  point  is  maximum  for  its  row 
(input  only  1st  and  2nd  points) 


NOTE:  In  specifying  a  circular  arc,  the  points  are  ordered 

such  that  a  counterclockwise  direction  about  the  center 
is  obtained  upon  moving  along  the  boundary. 
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DATA  NO.  5  -  BOUNDARY  CONDITION  CARDS 

Each  card  assigns  a  particular  boundary  condition  to  a  block  of 
elements  bounded  by  11 .  I2,  J1 ,  and  J2.  For  a  line,  II  =  12,  or 
J1  =  J2.  For  a  point,  II  =  12,  and  J1  =  J2. 

Format  (415,  110,  3F10.0) 

Columns  1-5  Minimum  I 

6-10  Maximum  I 

11-15  Minimum  J 

1 6-20  Maximum  J 

21-30  Boundary  condition  code 

31-40  Radial  boundary  condition  code,  XR 

41-50  Axial  boundary  condition,  XZ 

51-60  Tangential  boundary  condition,  XT 

If  the  number  in  Columns  21-30  is 

0  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-load,  and 
XT  is  the  specified  0-load. 

1  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-load,  and 

XT  is  the  specified  0-load. 

2  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  0-load. 
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3  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  e-load. 

4  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-load,  and 

XT  is  the  specified  e-displacem^^nt. 

5  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-load,  and 

XT  is  the  specified  e-displacement. 

6  XR  is  the  specified  R-load, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  e-displacement. 

7  XR  is  the  specified  R-displacement, 

XZ  is  the  specified  Z-displacement,  and 
XT  is  the  specified  e-displacement. 
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DATA  NO.  6  -  MATERIAL  BLOCK  ASSIGNMENT  CARD 


Each  card  assigns  a  material  definition  number  to  a  block  of 

elements  defined  by  the  I  and  J  coordinates. 

Format  (515,  2F1 0.0,  215) 

Columns  1-5  Material  definition  number 

(1  through  6) 

6-10  Minimum  1 

11-15  Maximum  I 

1 6-20  Minimum  J 

21-25  Maximum  J 

26-35  Material  principal  property 

inclination  angle  BETA  which  defines 
N-S  plane  orientation  relative  to 
z  direction  (see  Figure  4) 

36-45  Material  principal  property 

inclination  angle  ALPHA  which  defines 
the  orientation  of  the  N-T  plane 
relative  to  the  r-z  plane  (see  Figure  4) 

46-50  lANG  (If  lANG  =  0.  then  ALPHA  is 

same  for  total  material  block. 

If  lANG  =  1 ,  the  ALPHA  varies  in 
sign  in  the  I  direction  from 
element  to  element  every  NANG 
elements.  This  will  allow  for 
equal  but  opposite  helical  angles.) 

51-55  NANG  (Number  of  elements  in  the  I 
direction  with  the  same  ALPHA) 

56-60  IJDIR  (Gap  direction  0,1,  or  2) 

NOTE:  Gap  elements  will  be  identified  by  a  material  definition 
number  of  2. 
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DATA  NO.  7  -  PLOT  TITLE  CARD 


Forniui  (20A4) 

Columns  1-80  Title  (Title  printed  under  each  plot) 


DATA  NO.  8  -  PLOT  GENERATION  INFORMATION  CARD 
Format  (2F10.0) 

Columns  1-10  RMAX  (Maximum  r  coordinate  of  mesh) 

1 1-20  ZMAX  (Maximum  z  coordinate  of  mesh) 
NOTE;  Use  only  if  IPLOT  =  1  (plot  required). 
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DATA  NO.  9  -  TEMPERATURE  FIELD  INFORMATION  CARDS 


If  NUMTC  in  columns  6-10  of  the  CONTROL  CARD  is  greater  than 
1,  tlie  temperaiUic  ■'ielcJ  is  given  on  cards.  One  Card  must 
be  supplied  for  each  point  for  which  a  temperature  is 
specified. 

Format  (3F10.0) 


Columns 

1-10 

R  coordinate 

11-20 

Z  coordinate 

21-30 

Temperature 

If  NUMTC  in  columns  6-10  of  the  CONTROL  CARD  is  -2,  a  constant 
temperature  field  is  specified:  the  value  is  given  on  a 
single  card. 

Format  (FI  0.0) 

Columns  1-10  Temperature 
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DATA  NO.  1&  -  MATERIAL  PROPERTY  INFORMATION  CARDS 

The  following  group  of  cards  must  be  specified  for  each 
matertal  (niaxirr.uni  oi  S). 


MATERIAL  IDENTIFICATION  CARD 

Format  (215,  2F10.0) 

Columns  1-5 

Material  identification  number 

6-10 

Number  of  temperatures  for  which 
properties  are  given  (12  maximum) 

11-20 

Mass  density  of  material  (if 
required) 

21-30 

Thermal  expansion  parameter  (If 

1 ,  free  thermal  expansion  on  the 
material  property  cards;  othenwise, 
coefficients  of  thermal  expansion 
are  on  the  material  property  cards.) 

MATERIAL  PROPERTY  CARDS 

Two  cards  are  required  for  each  temperature. 

First  Card 

Format  (7F10.0) 

Columns  1-10 

Temperature 

11-20 

Modulus  of  elasticity,  E^ 

21-30 

Modulus  of  elasticity,  Eg 

31-40 

Modulus  of  elasticity,  Ej 

41-50 

Poisson’s  ratio,  v^g 

51-60 

Poisson’s  ratio, 

61-70 

Poisson’s  ratio,  Vg^ 
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Second  Card 


Format  (6F10.0) 
Columns  1-10 
11-20 
21-30 
31^0 
41-50 
51-60 


Shear  modulus,  G„js 
Shear  modulus,  Ggj 
Shear  modulus,  G-^n 
afj  or 
OgT  or  Os 
ajT  or  a-p 
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DATA  NO.  1 1  -  YIELD  STRESS  CARDS 


(not  needed  in  elastic  version) 

Format  (7F10.0) 

Columns  1-10  Yield  stress  in  tension  in  N 

direction 

1 1-20  Yield  stress  in  tension  in  S 
direction 

21-<30  Yield  stress  in  tension  in  T 
direction 

31-40  Yield  stress  in  shear  in  NS 
direction 

41-50  Yield  stress  in  shear  in  NT 
direction 

51-60  Yield  stress  in  shear  in  TS 

direction 

61-70  Hardening  parameter;  C 
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DATA  NO.  12  -  INERTIA  LOAD  CARD 


Format  (3F10.0) 

Starting  with  this  input  card  and  including  the  boundary 
force  cards,  this  data  is  to  be  input  as  a  block  for 
each  load  step,  that  is,  NLINC  times.  There  are  the 
following  exceptions  to  this: 

(a)  If  INERT  =  0,  then  this  card  is  to  be 
omitted  completely  (no  inertia  load). 

(b)  If  INCI  =  0,  then  this  card  is  not 
repeated,  but  appears  in  first  block 
only  (the  inertia  loads  are  constant 
for  each  load  step). 

(c)  If  INCF  =  0,  then  the  following  boundary 
pressure  and  shear  cards  are  to  be  given 
only  for  the  first  block  and  not  repeated 
again  (the  pressure  and  shear  loads  are 
constants  for  each  load  inaement). 


Columns  1-10 

ACELZ  (axial  acceleration) 

11-20 

ANGVEL  (angular  velocity) 

21-30 

ANGACC  (angular  acceleration) 
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DATA  NO.  13  -  BOUNDARY  PRESSURE  CARDS 

One  card  is  required  for  each  boundary  element  which  is 
subjected  to  a  normal  pressure;  that  is.  the  number  of 
these  cards  is  NUMPC  for  each  load  increment. 

Format  (215,  F10.0) 

Columns  1-5  Nodal  point  M 

6-1 0  Nodal  point  N 

1 1-20  Normal  pressure 

As  shown  in  Figure  A-1,  the  boundary  element 
must  be  on  the  left  when  progressing  from  M  to  N. 
Surface  normal  tension  is  input  as  a  negative  pressure. 


Figure  A-1 .  Example  of  Boundary  Pressure  Loading 
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DATA  NO.  14  -  BOUNDARY  SHEAR  CARDS 


One  card  is  required  for  each  boundary  element  which  is 
subjected  to  surface  shear;  that  is,  the  number  of 
these  cards  is  NUMSC  for  each  load  increment. 


Format  (215,  F10.0) 

Columns  1-5 

Nodal  point  M 

6-10 

Nodal  point  N 

11-20 

Surface  shear 

As  shown  in  Figure  A-2.  the  boundary  element 
must  be  on  the  left  when  progressing  from  M  to  N. 
The  positive  sense  of  the  shear  is  from  M  to  N. 


Figure  A-2.  Example  of  Boundary  Shear  Loading. 
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DATA  NO.  15  -  BOUNDARY  TRANSVERSE  SHEAR  CARDS 

One  card  is  required  for  each  boundary  element  which  is 
subject  to  the  transverse  shear;  that  is,  the  number  of  these 
cards  is  NUMSC  for  each  load  increment. 


Format  (2I5,  FI 0.0) 

Columns  1-5 

Nodal  point  M 

6-10 

Nodal  point  N 

11-20 

Surface  transverse  shear 

Figure  A-3.  Example  of  Boundary  Transverse  Shear  Loading. 
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APPENDIX  B: 

COMPUTER  LISTING  OF  SUBROUTINE  ELPLSS 
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SUBROUTINE  ELPLSS(MTYPE) 

COMMON/ 1 NCR/NOL , I NER  T , NUMMAT , S I GTOT ( 12 , 1 00 ) 
l,EPSTOTn2,100) 

C0MM0N/ARG1/SIC1(18) .EPS1(18) ,DEPSP(12) .CEPSP(6,6) 
COMMON/PLAS/ALFA( 6 , 100 ) ,S I GYLO( 7 ,6 ) , I FGPL { 100 ) 
C0MM0N/ARG/RRR(5) ,ZZZ(5) ,RR(U) ,Z2(4) ,S( 15,15) ,P{ 15i .TT(6) , 
1H(6,15) .CRZ(6,6) ,XI (10) .ANGLE(4),SIG(18) ,EPS(18),N 
C0MM0N/RESULT/BS(6,15),D(6,6),C(6,6),AR,BB(6,9),CNS(6,6) 
DIMENSION  DFDSIG(6) ,SIGYB{3) ,DUHH(6),D1(6,6) 

C  FIND  DF/DSIG 

S I GYB ( 1 ) =1 . /S I GYLD ( 1 ,MTYPE )**2- 1 ./S I CYLD (2 ,MTYPE )**2 
1  -1 ./SIGYLD(3,MTYPE)**2 

S I G YB ( 2 )  =  1 . /S I GYLD ( 2 , MT YP  E ) **2- 1 . /S I GYLD ( 1 , MTYP  E ) **2 
1  -l./SIGYLD{3,HTYPE)**2 

S  I GYB  ( 3 )  =  1 . /S  I  GYLD  ( 3  ,  MTYPE )  ••*2- 1 . /S I  GYLD  ( 2 ,  KTYPE )  **2 
1  -l./SIGYLD(l,MTYPE)**2 

DO  105  1=1,6 

SIG1( I )=SIGTOT( l+6,N)-ALFA( I ,N) 

105  CONTINUE 

DFDSIG(1)=2.*SIG1(1)/S1GYLD(1,MTYPE)**2+SIGYB(2)*SIG1(3) 

1  +SIGYB(3)*SIG1(2) 

DFDSIG(2}=2.*SIG1(2)/SIGYLD{2,MTYPE)**2  +5 IGYB{ 1 )*S IG1 ( 3 ) 

1  +SIGYB(3)*SIG1{1) 

DFDSIG(3)=2.*SIG1(3)/SIGYLD(3,MTYPE)**2+SIGYB(1)*SIG1(2) 

1  +SIGYB{2)*SIG1( 1) 

DO  110  1=4,6 

110  DFDSIG( I )=2.*SIG1( I )/SIGYLD( I ,MTYPE)**2 
DO  120  1=1,6 
DUMM( I )=0.0 
DO  120  J=1,6 

120  DUMM( I )=DUMM( I )+CNS( I , J)*DFDS IG( J) 

D=0. 

DO  125  1=1,6 

125  D=D+DFDSIC( I )*DUMM( I ) 

DO  126  1=1,6 

126  D=D+SIGYLD(7,KTYPE)*DFDSIG( l)**2 
D=1,/D 

DO  130  1=1,6 
DO  130  J=1,6 

130  D1( I ,J)=DFDSIG( I )*DFDSIG( J) 

DO  140  1=1,6 
DO  140  J=1,6 
CEPSP( J, I )=0.0 
DO  140  K=1,6 

140  CEPSP( J,l)=CEPSP( J, I )+D*CNS( I ,K)*D1(K,J) 

DO  150  1=1,6 
DO  150  J=1,6 
D1( I ,J)=0.0 
DO  150  K=1,6 

150  D1{I,J)  =D1(I,J)‘*'CNS(I,K)*CEPSP{K,J) 

DO  160  1=1,6 
DO  160  J=1,6 

160  CNS( I ,J)=CNS( I ,J)-D1{ I ,  J) 

RETURN 

END 
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2  Administrator  1 

Defense  Technical  Info  Center 
ATTN;  DTIC-DDA 
Cameron  Station 
Alexandria,  VA  22304-6145 

1 

1  Commander 

U.S.  Army  Materiel  Command 
ATTN:  AMCDRA-ST 
5001  Eisenhower  Avenue 
Alexandria,  VA  22333-0001 

1 

1  Commander 

U.S.  Army  Laboratory  Command 
ATTN:  AMSLC-DL 
2800  Powder  Mill  Road 

Adelphi,  MD  20783-1145  1 

2  Commander 

U.S.  Army  Armament  Research, 

Development,  and  Engineering  Center 
ATTN;  SMCAP.  !MI1  taa««.oniy)i 

Picatinny  Arsenal,  NJ  07806-5000 

2  Commander 

U.S.  Army  Armament  Research, 

Development,  and  Engineering  Center  (Unciaas.  oniy)i 
ATTN:  SMCAR-TDC 
Picatinny  Arsenal,  NJ  07806-5000 

1  Director 

Benet  Weapons  Laboratory  1 

U.S.  Army  Armament  Research, 

Development,  and  Engineering  Center 
ATTN:  SMCAR-CCB-TL 
Watervliet,  NY  12189-4050 


Commander 

U.S.  Army  Missile  Command 
ATTN;  AMSMI-RD-CS-R  (DOC) 

Redstone  Arsenal,  AL  35898-5010 

Commander 

U.S.  Army  Tank- Automotive  Command 
ATTN;  ASQNC-TAC-DIT  (Technical 
Information  Center) 

Warren,  Ml  48397-5000 

Director 

U.S.  Army  TRADOC  Analysis  Command 
ATTN;  ATRC-WSR 

White  Sands  Missile  Range,  NM  88002-5502 
Commandant 

U.S.  Army  Field  Artillery  School 

ATTN;  ATSF-CSI 

Ft.  Sill,  OK  73503-5000 

Commandant 

U.S.  Army  Infantry  School 

ATTN:  ATSH-CD  (Security  Mgr.) 

Fort  Benning,  GA  31905-5660 

Commandant 
U.S.  Army  Infantry  School 
ATTN;  ATSH-CD-CSO-OR 
Fort  Benning,  GA  31905-5660 

Air  Force  Armament  Laboratory 

ATTN;  WUMNOI 

Eglin  AFB,  FL  32542-5000 

Aberdeen  Proving  Ground 


(Unclass.  oniy)i  Commander 

U.S.  Army  Armament,  Munitions 
and  Chemical  Command 
ATTN:  AMSMC-IMF-L 
Rock  Island,  IL  6t  299-5000 


2  Dir,  USAMSAA 
ATTN:  AMXSY-D 

AMXSY-MP,  H.  Cohen 

1  Cdr,  USATECOM 
ATTN:  AMSTE-TC 


1 


Director 

U.S.  Army  Aviation  Research 
and  Technology  Activity 
ATTN:  SAVRT-R  (Library) 
M/S  219-3 

Ames  Research  Center 
Moffett  Field,  CA  94035-1000 


3  Cdr.  CRDEC,  AMCCOM 
ATTN:  SMCCR-RSP-A 
SMCCR-MU 
SMCCR-MSI 

1  Dir,  VLAMO 

ATTN:  AMSLC-VL-D 


10  Dir,  BRL 

ATTN:  SLCBR-DD-T 
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2  Director 
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4  Commander 

U.S.  Army  Armament  Research, 
Development,  and  Engineering  Center 
ATTN;  SMCAR-CCH-T, 

S.  Musalli 
SMCAR-CC, 

J.  Hedderich 
E.  Fennell 
R.  Price 

Picatinny  Arsenal,  NJ  07806-5000 

1  Commander 

U.S.  Army  Armament  Research, 
Development,  and  Engineering  Center 
ATTN:  SMCAR-TD,  T.  Davidson 
Picatinny  Arsenal,  NJ  07806-5000 

2  Commander 

U.S.  Army  Materials  Technology  Laboratory 
ATTN:  SLCMT-MEC, 

B.  Halpin 

T.  Chou 

Watertown,  MA  02172-0001 
2  Director 

Lawrence  Livermore  National  Laboratory 
ATTN:  S.  DeTeresa 

R.  M.  Christensen 
P.O.  Box  808 
Livermore,  CA  94550 

4  Director 

Sandia  National  Laboratories 
Applied  Mechanics  Department, 
Division-8241 
ATTN:  C.  W.  Robinson 
G.  A.  Benedetti 
K.  Perano 
W.  Kawahara 
P.O.  Box  969 

Livermore,  CA  94550-0096 
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2  Battelle  Pacific  Northwest  Laboratory 
ATTN:  M.  Smith 
M.  Garnich 
P.O.  Box  999 
Richland,  WA  99352 

1  Director 

Los  Alamos  National  Laboratory 

ATTN:  D.  Rabern 

WX-4  Division,  Mail  Stop  G-787 

P.O.  Box  1663 

Los  Alamos,  NM  87545 

2  David  Taylor  Research  Center 
ATTN:  R.  Rockwell 

W.  Phyillaier 

Bethesda,  MD  20054-5000 

2  Olin  Corporation 

Flinchbaugh  Operations 
200  E.  High  Street 
ATTN:  B.  Stewart 
E.  Steiner 

Red  Lion,  PA  17356 

2  Alliant  Tech  Systems 
5640  Smetana  Drive 
ATTN:  J.  Bode 
K.  Ward 

Minnetonka,  MN  55343 

2  Project  Manager 

Tank  Main  Armament  Systems 
ATTN:  SFAE-AR-TMA-MD, 

C.  Kimker 
SFAE-AR-TMA-ME, 

K.  Russell 

Picatinny  Arsenal,  NJ  07806-5001 

1  Zak  Technologies,  Inc. 

ATTN:  Dr.  Adam  R.  Zak 
2310  Belmore  Dr. 

Champaign,  IL  61821 


54 


